import sys, pylab
sys.path.append( '../lib/' )
import numpy as np
from gtd import *

# Fixed Variables
kB  = 1.3806488e-23
c = 3e8
hw = 13.152
h0 = 10.0
R0 = 13.6
nu0 = 148e9
dnu = 20e9
T = 300

# Derived variables
lambd = c/nu0
phi  = 2*pi - np.arctan( R0/(hw-h0) )
s    = np.sqrt( R0**2 + (hw-h0)**2 )
k    = 2*pi/lambd

def Bbflux (T, nu, dnu):
    """
    @brief Blackbody flux
    """
    const=(nu/c)**2
    B = 2*kB*T*const*dnu	
    return B

def Fwall (phi_p):
    """
    @brief Evaluates the incident flux density at a given incidence
           angle "phi_p"
    """
    B = Bbflux(T, nu0, dnu)
    F = (R0/hw)*np.cos(phi_p) + np.sin(phi_p)
    return B*F


#Integration variables
phi_p_min=0.0
phi_p_max=pi/2
n_int=1000


vector_phi_p = np.linspace(phi_p_min,phi_p_max,n_int)
delta=(phi_p_max-phi_p_min)/n_int
vector_Fwall = Fwall(vector_phi_p)

# Obtain the flux density contribution from every incident angle at the telescope
# after difracting over the ground shield.
# Note that EdiffPlane acts over the electric field and "vector_Fwall" contains the
# a measure of the power proportional to E^2, so we need to provide sqrt(vector_Fwall)
# "vector_D" contains the electric field intensity at the telescope.
vector_D = np.abs(EdiffPlane( np.sqrt(vector_Fwall), s, k, phi, vector_phi_p ))

# Here we get the integral by summing all the constributions times the fractional angle
# delta. The 2*pi factor comes from integrating in azimuth for all directions
F_int = 2*pi*sum(vector_D**2)*delta 

# We can integrate that flux through the aperture of the telescope and attenuate it by
# sidelobe power (~-80 dB) to estimate the loading effect on the detectors
aperture = 19.3 # m^2
sl_level = [-79.7, -80.4, -81.3, -84.9, -87.6, -89, -88.6] # dB
elev = [31, 35, 40, 45, 50, 55, 60]
power = F_int*aperture*np.power(10,np.array(sl_level)/10)

print "Flux/A = %7.3g [W/m^2]"%F_int
for e in xrange(len(elev)):
    print "Loading (el = %d) = %7.3e [pW]"%(elev[e],power[e]*1e12)

pylab.plot(elev, power*1e15, "o:")
pylab.xlabel("Elevation [deg]")
pylab.ylabel("Ground Loading [fW]")
pylab.savefig("ground_loading_vs_el.png")
